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A  TRIANGULAR  THIN  SHELL  ELEMENT 
FOR  THE  LINEAR  ANALYSIS  OF 
STIFFENED  COMPOSITE  STRUCTURES 

By 
Edward  Wilson  and  Gilles  Cantin 

ABSTRACT 

A   three  node  flat  shell  element  with  six  engineering  displacement 
degrees-of-f reedom   at   each   node   is   developed.     The   basic 
formulation   allows   for   the   arbitrary  location  of  the  reference 
surface  in  which  the  membrane  forces  and  bending  moments  are  fully 
coupled. 

The  well-known,  highly  accurate,  DKT  bending  element  is  combined 
with  a  higher  order  membrane  element  in  order  to  obtain  a 
consistent  formulation.  The  higher  order  membrane  behavior  is 
obtained  by  the  introduction  of  three  additional  normal  rotational 
degrees-of-f reedom. 

This  report  presents  a  summary  of  the  theoretical  steps  involved 
in  the  development  of  the  element.  The  accuracy  of  the  element  is 
illustrated  by  the  solution  of  several  standard  problems  and  a 
comparison  of  results  with  other  thin  shell  elements.  The  FORTRAN 
77  listing  of  the  subroutines  which  form  the  basic  element 
matrices  contain  less  than  300  statements  and  is  presented  in 
order  to  illustrate  that  the  computer  implementation  of  the 
element   is   relatively  simple. 


-  11  - 


TABLE  OF  CONTENTS 

ABSTRACT  i 

TABLE  OF  CONTENTS  ii 

INTRODUCTION  1 

BASIC  EQUATIONS  3 

BENDING  APPROXIMATION  -  THE  DKT  ELEMENT  6 

MEMBRANE  APPROXIMATION  10 

COMPUTER  PROGRAM  IMPLEMENTATION  12 

NUMERICAL  EXAMPLES  12 

Cantileaver  Beam  Example  13 

Twisted  Beam  14 

Scordelis-Lo  Roof  15 

Spherical  Shell  16 

FINAL  REMARKS  17 

REFERENCES  17 

FORTRAN  LISTING  OF  18  DOF  TRIANGULAR  SHELL  ELEMENT  18 


-  1  - 


INTRODUCTION 

General  Background 

The  use  of  composite  materials  allows  for  the  efficient  design  of 
many  different  types  of  structural  systems.  One  of  the  major 
advantages  of  the  material  is  that  different  stiffness  and 
strength  properties  can  be  obtained  in  different  directions. 
Therefore,  more  efficient  structures  can  be  obtained  since  the 
material  can  be  concentrated  in  the  directions  of  maximum 
stresses. 

Most  of  the  existing  finite  element  programs  do  not  have 
sufficient  generality  to  consider  such  material  properties.  Also, 
in  the  case  of  thin  shell  structures  very  few  programs  have  the 
ability  to  consider  shells  in  which  the  bending  and  membrane 
forces  are  coupled.  In  addition,  problems  associated  with  the 
modelling  of  complex  shell  structures  with  thin  shell  elements 
exist  since  the  classical  thin  shell  formulation  does  not  have 
stiffness  terms  associated  with  the  normal  rotational  degrees  of 
freedom.  Therefore,  the  user  of  the  program  is  often  required  to 
add  artificial  members  to  a  finite  element  model  in  order  to  avoid 
numerical  instability  in  the  solution  of  the  finite  element 
system.  The  purpose  of  this  report  is  to  present  a  new  thin  shell 
element  which  is  sufficiently  robust  to  solve  the  above  mentioned 
problems . 

Recent  Research 

Within  the  past  two  years  several  papers  have  presented  methods 
which  introduce  a  normal  rotation  in  order  to  improve  the  membrane 
behavior  of  plane  elements.  Carpenter,  Stolarski  and  Belytschko 
present  a  flat  triangular  shell  element  with  improved  membrane 
interpolation.  [1]  They  introduced  normal  mid-side  displacements 
of  the  constant  strain  triangle.  The  normal  displacements  are 
eliminated  and  node  rotations  are  introduced  by  the  use  of  a  cubic 
constraint  function  along  each  side  of  the  triangle.  A  one  point 
integration  method  is  used  in  order  to  eliminate  membrane  locking 
within  the  elements.  The  element  yields  very  accurate 
displacements;  however,  the  element  is  rank  deficient  and  is 
unstable   for   certain  geometries. 
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Taylor  and  Simo  applied  the  same  basic  approach  as  presented  in 
reference  [1]  to  improve  the  membrane  behavior  of  quadrilateral 
elements  [2] .  For  many  problems  excellent  displacements  and 
stresses  are  obtained.  However,  for  shell  structures  such  as  a 
twisted  beam  the  displacements  become  very  large  as  the  mesh  was 
refined.  In  addition,  the  flat  quadrilateral  element  does  not 
accurately  model  many  common  types  of  shell  geometries.  Also,  the 
DKQ  formulation  was  used  to  form  the  bending  stiffness  which  has 
proven   to   be   not   as   accurate   as   the   DKT   formulation. 

Bergan  and  Felippa  have  developed  a  triangular  membrane  element 
with  normal  rotational  degrees-of-f reedom  [3] .  The  formulation  is 
based  on  the  "free  formulation"  and  uses  the  continuum-mechanics 
definition  of  rotation.  The  element  passes  the  patch  test  and  is 
stable  for  all  applications.  The  element  produces  good  values  of 
displacements;  however,  the  values  for  stresses  are  poor  compared 
to   the   values   obtained   from   the   Taylor   quadrilateral   [4] . 

The  three  elements  previously  mentioned  have  not  been  used  for 
thin  shells  in  which  the  membrane  and  bending  forces  are  coupled. 
Therefore,  one  of  the  purposes  of  this  report  is  to  develop  an 
element  which  has  a  consistent  formulation  for  both  the  bending 
and  membrane  behavior.  In  addition,  the  problems  with  the 
instability  associated  with  normal  rotational  degree-of-f reedom 
will  be  studied  and  a  simple  technique  is  suggested  in  order  to 
avoid   this   problem. 
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BASIC  EQUATIONS  -  ORTHOTROPIC  MATERIALS 

The  18  x  18  triangular  shell  element  stiffness  matrix  for  a  stiffened 
composite  material  as  shown  in  figure  1  can  be  directly  calculated 
from   the   following   well-known   equation: 


K  = 


BT  DBdA 


(1) 


The  6x6  D  matrix  relates  the  forces  to  the  deformations  which  are 
associated  with  a  differential  element  of  area  dA.  Including  thermal 
deformations  the  force-deformation  relationship  can  be  expressed  by 
the   following   matrix   equation: 


f  =  D  £   +   fo 


where 


and 


f_T  =  [  mi  i  nit  mi  x  f  i  i  fix  fix  ] 


e_T  =  [  ki  i  kxx  kix  tn  cxx  eix  ] 


(2) 


(3) 


(4) 


The   positive   definition   of   these   forces   and   deformations   is 
illustrated   in   figure   2. 


Normally  the  matrix  D  cannot  be  defined  directly  for  complex 
materials.  However,  the  inverse  D~ l  can  normally  be  easily  calculated 
from  the  basic  principles  of  mechanics  or  determined  experimentally. 
Therefore,  the  terms  for  D-1  are  normally  specified  as  input  to  a 
computer  program.  The  numerical  values  of  D  are  then  evaluated  within 
the  element  stiffness  subroutine.  Hence,  the  basic  behavior  of  the 
thin  shell,  including  thermal  deformations,  is  expressed  in  the 
following   form: 
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Figure  1.  EXAMPLE  OF  ANISOTROPIC  SHELL 
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Figure  2.  DEFINITION  OF  POSITIVE  FORCES 
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Or  in  terms  of  matrix  notation: 

L  ~   D"1  f   +   cp  (5) 

Where  dT  is  the  temperature  change  and  cu  to  a*  are  the  measured 
thermal  expansion  coefficients.  Hence,  the  thermal  forces  indicated 
in   equation   (2)   are   calculated   from: 

f o  =  -  D  e.o  (6) 

Each  flexibility  term  in  equation  (6)  has  a  direct  physical  meaning. 
For   example,   the  term  Ci a  is  the  curvature  k\i   due  to  a  unit  force, 
ftz      =   1.   Whereas,   the  term  C2 1  is  the  strain  tn    at  the  reference 
plane   caused  by  the  application  of  a  unit  bending  moment,  mi  1  =1. 
It  is  apparent  that  these  terms  can  be  best  determined  experimentally 
for   complex  composite  materials.  Also,  the  values  of  these  terms  are 
dependent   on   the   definition  of   the  reference  plane  which  must  be 
defined   at   the   same  time  as  the  flexibility  terms  are  determined. 
For   the   special   case   of   constant   thickness  isotropic  shells  the 
mid-surface   is   the   logical   reference   plane  and  the  terms  Cij  and 
several   other   flexibility   terms   are   zero   in   equation   (5) . 

In  equation  (1)  the  6x18  B  matrix  defines  the  relationship  between 
the  deformation  terms  and  the  node  displacements  v  in  the  local 
1,2,3   system.   Or,   in  matrix   form: 

£_   =   B   v  (7) 

which   can   be   written  in   submatrix   form   as 

1p   =   Bj»   v  (8a) 

£■   =   Ba   v  (8b) 

where  the  "p"  and  "mH  indicate  the  plate-bending  and  membrane  terms 
respectively. 
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BENDING  APPROXIMATION  -  THE  DKT  ELEMENT 

The  development  of  the  BP  matrix  is  based  on  the  standard  DKT  elemen 
[5] .  Because  the  bending  and  membrane  behavior  are  coupled  the  DK 
formulation  will  be  summarized  here  in  order  that  the  Ba  matrix  wil 
be  developed  with  consistent  approximations  in  the  next  section  o 
this   report. 

The  DKT  element  is  based  on  the  independent  expansion  of  the  inplan 
rotations  of  the  reference  surface  for  a  6  node  triangle  which  i 
shown  in  figure  3.  If  the  local  1  and  2  directions  are  indicated  b 
the  local  x  and  y  coordinates  the  finite  element  approximation  i 
written   in   the   following   form: 

Bz    =   ai    +   at  x  +   as  y   +      fca*  x2    +      as  xy   +   J4a*  y1 

(9 
Br    =   a7    +   as x  +   a* y  +  Jiai  ox1    +  anxy  +  fcai  s y* 

The      six      constants   ai    to   ac ,    ax ,    can  be   expressed   in   terms   of    the   si 
node      rotations      $x i    to   Bz»    by   an   inversion   of    a    6x6   matrix  which  wi! 
produce      an      equation      of      the      following      form: 

ax    =   H   B_x  (10) 

The  same  6x6  matrix,  H,  will  relate  the  constants  a?  to  an,  ar  ,  t 
the   node   rotations   B_t  . 

From   the   theory  of   thin  plates   the   curvature-displacement 
relationships   are  defined  by   the   following  equations: 

kz  i    =   Bi  , .      =     ax    +  a*  x  +  as  y  ( 

k*F    =  By  ,w      =     ■•    +  anx  +  aixy  ( 

kz  f    ■  3z  , t    +   (J*  ,x 
.a 

-at+asx  +  aty  +  at+aiox+any  ( 

Or,  in  the  following  matrix  form: 

tf    =  Gp  a  (12! 


-  7  - 


^►X 


(a)  BASIC  ROTATIONAL  UNKNOWNS 


Ljj  is  the  length  from  I  to  J 


(b)  DISPLACEMENTS  ALONG  TYPICAL  SIDE 


Figure  3.   DISPLACEMENT  APPROXIMATION  -  DKT  ELEMENT 
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Where 


apT    =    [a*    +   a3    +   as    +   a«    +   as    +   a*    + 
as    +   a»    +   ai  o    +   an    +   an    ] 


and 


Gp  = 


lOxyOOOOOO 
OOOOOOlOxy 
OlOxylOxyO 


(13 


(14 


The 


six  rotational  degrees-of-f reedom  which  are  associated  withi 
typical  side  I-J  of  the  element  are  shown  in  figure  (3a) .  The< 
rotations  can  be  transformed  to  a  n-s  reference  system  which  { 
parallel  and  normal  to  a  typical  element  side  as  shown  in  figure  (3b, 
The   basic   DKT   constraints   are   enforced   as   follows: 


1.   The  mid-side  rotation  9.  .  is  set  to  the  average  of  the  valuj 
at  point  I  and  J.   Or, 


68.  =  (  6»i  +  e.j  )  /  2 


(15 


2.   The  s-displacements,  above  and  below  the  reference  plane,  ac 
the  normal  displacement  in  the  z-direction  are  forced  to  be 
cubic  functions  since  the  transverse  shear  strain  is  in  the 
s-direction  is  set  to  zero.   Therefore,  the  mid-side  normal 
rotation  must  satisfy  the  following  equation: 


enxa  -   -  (e.i  +6.j)/4   -   3(wi  +wj)/(2Lu) 


(16 


The  constraint  specified  by  equation  (15)  will  force  the  norm! 
displacement  along  the  element  sides,  above  and  below  the  referenE 
plane,  to  be  linear  function.  Hence,  displacement  and  slope 
compatibility  is  satisfied  along  the  sides  of  all  elements.  Since  : 
attempt  is  made  to  set  the  transverse  shear  strains  to  zero  within  t: 
element  the  name  Discrete  Kichoff  Triangle  was  selected  as  the  name  I 
this   element. 
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Equation  (15)  and  (16)  can  be  summarized  in  matrix  form  as 


8s 

1/2         0 

6s  i 

1/2 

0 

9s  j 

9* 

0      -1/4 

On  i 

0 

-1/4 

Oh  j 

+  1.5/L 


i  j 


0 
-1 


J_ 


Wi 
Wj 


The  relationship  between  the  N-S  and  X-Y  coordinate  systems  are 


8s 

eN 

and 
9x 
9y 


Cos5  SinS 

-Sin5  Cos5 

Cos5  -Sin5 
SinS    Cos5 


9x 
9y 


6s 
On 


Equation  (17)  can  now  be  written  in  the  X-Y  system  as 


9x 

9y 

Where 


9xi 
9yi 


Tl 

T2 

9xj 

+    1.5/Lu 

-Ts 

Ts 

Wi 

T2 

T3 

9y  j 

Tc 

-Tc 

Wj 

(17) 


(18a) 


(18b) 


(19) 


Tl  =  0.5  Cos2 5  -  0.25  Sin2 6 

T2  ■  1.5  Sin5  Cos5 

T3  =  0.5  Sin2 5  -  0.25  Cos2 6 

Ts  =  1.5  Lu  Sin5 

Tc  =  1.5  Li  j  Cos5 

The  six  rotational  degrees-of-freedom  associated  with  the  three 
mid-side  nodes  of  the  triangle  can  be  eliminated  by  the  application  of 
equation  (19) .  These  transformations  can  be  summarized  by  a  matrix 
equation   of   the   following   form: 


e  =  tp  w 


(20) 


where   Tp   is   a  12x9  matrix  and  W  represents  a  9x1  vector  of  9x  i  ,  9v  i 
and  wi   for   the   three   nodes   of   the   triangle. 
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MEMBRANE  APPROXIMATION 

The   membrane  behavior  of  the  triangular  shell  element  is  based  on  th 

basic   six  node  quadratic  membrane.   If  the  12  coefficents  are  define 

by   the   12   x  1  vector  b  the  membrane  strains  can  be  expressed  in  th 
following   form: 


£h  =  Gh  b 


(21) 


As  in  the  case  of  the  DKT  plate  bending  element  the  three  mid-sid 
displacements  are  rotated  to  the  local  N-S  coordinate  of  each  side  a 
shown  in  figure  (4).  In  order  to  maintain  displacement  compatibilit 
between  element  the  displacement  us  is  assumed  to  be  linear  along  eac 
side  and  the  displacement  un  is  a  cubic  function.  These  assumption 
can  be  summarized  by  the  following  equations  for  the  displacements  a 
the   mid-side   nodes: 


us  =  ( us  i  +  us  j  )  /  2 

un    =    (uni     +    unj)/2      +   Lu(6zi     -   ©zj)/8 


(22) 


These  equations  can  be  written  in  terms  of  the  X-Y  coordinate  system 
as 


ux 
uy 


.5 

0 


o" 

Ux  I 

+ 

".5 

o" 

Tux  j 

5 

Uy  i 

0 

.5 

1  Uy  j 

125Lu 


-SinS   Sin5 
Cos5  -Cos5 


e2i 

e2J 


(23) 


The  six  translational  displacements  at  the  midside  nodes  can  now  b 
eliminated  and  three  rotational  unknowns  are  introduced  at  th 
vertices  by  the  direct  application  of  equation  (23)  .  The  same  basi 
approach  which  was  used  in  the  DKT  formulation  is  now  applied  in  orde 
to   form   the  matrix  equation  of   the   following   form: 


b  =  Tm  u 


(24) 


Therefore,   the   three  membrane  strains  can  be  written  in  terms  of  th< 
nine   node   displacements   as 


Cm 


=   Gm   Tm 


U 


(25) 


It   is  now  possible  to  evaluate  the   24  x  24   element  stiffness  by  th< 
direct   application   of   equation   (1) . 
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Figure  5   TRIANGULAR  18  DEGREE-OF-FREEDOM  SHELL  ELEMENT 
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COMPUTER  PROGRAM  IMPLEMENTATION 

The  18  x  18  triangular  element  stiffness  matrix,  for  the  element  showi 
in   figure   5,   is   given   by   equation   (1) .       Where   the   the 
strain-displacement   matrix   B   can   now  be   written   in  terms  of  th 
bending   and  membrane   submatrices,   or 

Gp   0     Tp   0 
B  =  =   G(x,y)  T  (2 

0    Gm    0    T* 

Since   the   T  matrix   is   not  a  function  of  x  and  y  it  is  possible  t> 
rewrite   equation   (1)   in   the   following   form: 


K  ■  3?    G 


GT  D  G   dA   T  (2' 

The  matrix  G  is  very  sparse  and  contains  only  the  terms  1,  x  and  y 
therefore  the  integral  cab  be  evaluated  directly  with  a  minimum  o 
numerical  effort  as  illustrated  by  the  FORTRAN  listing  given  i: 
Appedix  A.  In  addition,  an  integral  reduction  factor  can  be  used  o: 
selective  terms  in  order  to  improve  the  membrane  performance  a: 
suggested   in  reference   (3) . 


NUMERICAL  EXAMPLES 

In  order  to  illustrate  the  behavior  of  the  element  and  to  compare  thi 
results   with  other  shell  elements  several  examples  will  be  presented 

Cantilever   Beam 

The  beam  shown  in  figure  6  is  idealized  by  a  1x4  rectangular  mesh  an< 
is  subjected  to  a  load  of  40  kips  at  the  tip  of  the  cantilever.  Th< 
theoretical  displacement  at  the  tip,  including  shearing  deformation 
is   0.3558   inches. 
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The  Taylor  quadrilateral  element  shell  yields  a  displacement  of  0.3467 
inches;  or  an  error  of  -1.02  percent.  Note  that  the  rotations  at  the 
base  of  the  cantilever  are  set  to  zero  which  is  inconsistent  with  the 
existence   of   shearing   deformations. 

The  element  presented  in  this  report,  TSHELL,  was  used  to  model  this 
beam  with  two  triangles  to  form  each  quadrilateral.  The  completely 
integrated  element  produces  a  displacement  of  0.2695  inches,  or  an 
error  of  -24.3  percent.  With  a  reduced  integration  factor  of  0.5  the 
tip  displacement   is   0.3726,   or   an   error   of   +4.5   percent. 

For   all   example   problems   presented   in   this   report   a   reduced 
integration  factor  of  0.5  is  used.   The  use  of  the  reduced  integration 
factor   has   the   major  advantage  over  one  point  integration  is  that  a 
"rank   deficiency"   is   not   introduced   into   the   element. 


Figure  6.   CANTILEVER  BEAM  EXAMPLE 
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Twisted  Beam 

The  twisted  beam  shown  in  figure  7  has  become  a  standard  test  problem 
for  thin  shell  elements  [6] .   Table  1  summarizes  the  results  obtained 
using  four  different  elements.   It  is  important  to  note  that  the 
SHELL  element  does  not  converge  as  the  mesh  is  refined.   The  reason 
for  this  unacceptable  behavior  is  because  the  element  is  flat  and  it 
cannot  model  the  twisted  surface  accurately.   The  new  triangular 
element  gives  good  results  for  this  problem. 


I2x  1.0 


WIDTH  =  I.I 
DEPTH  =  0.32 
TWIST  =  90° 
E  =  29.0x  I06 
POISSON  RATIO  =  0.22 
MESH  =  12  x2 

Unit  loads  at  end 


Figure  7.  TWISTED  BEAM  EXAMPLE 


Table  1.   Deflection  Under  Load  For  Twisted  Beam 


Element  Type  and  Mesh 

Reference  Values 

NASTRAN  QUAD4 
NASTRAN  QUAD 8 

SHELL  (2x12) 
(4x24) 

TSHELL  (2x12) 
(4x24) 


IN-PLANE 

error 

0.005424 

0 

0.005386 

-0.7 

0.005413 

-0.2 

0.005779 

+  6.4 

0.006849 

+26.2 

0.005390 

-0.6 

0.005399 

-0.5 

OUT-OF-PLANE 

erro 

0.001754 

0 

0.001727 

-1.5 

0.001750 

-0.2 

0.001993 

+13.7 

0.002750 

+56.9 

0.001717 

-2.1 

0.001735 

-1.1 
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Scordelis-Lo  Roof 

The  shell  structure  shown  in  figure  8  is  a  standard  test  problem  [6] . 
Table  2  summarizes  the  maximum  displacement  obtained  using  different 
elements  and  meshes.  For  this  problem  the  SHELL  element  results  are 
very  good.  However,  the  TSHELL  results  appear  to  converge  very 
slowly.  However,  for  even  the  coarse  mesh,  the  results  are  of 
acceptable   accuracy  for  normal   engineering  analysis. 


THICKNESS  =0.25 
E  =4.32x  I08 
POISSON  RATIO  =0.0 
LOADING:   90/unit  ores  in  z  dir. 
u  s  w  =  0  on  curved  edge 
MESH  NxN  on  quadrant 


Figure  8.   SCORDELIS-LO   CYLINDRICAL  SHELL 


Table  2.   Maximum  Displacement  Of  Scordelis-Lo  Roof 


Element  Type  and  Mesh 
Reference  Value 


SHELL 

N*4 

N=6 

N=8 

N=10 

TSHELL 

N=4 

N=8 

VALUE 
0.3086  ft. 


0.3036 
0.2982 


ERROR 

0.0  % 

+5.2  % 
+1.7  % 
+0.6  % 
+0.2  % 

-1.6  % 
-3.4  % 


N=12 


0.2997 


-2.9  % 
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Spherical  Shell 

A  spherical  shell  subjected  to  point  loads  is  shown  in  figure  9  and  is 
also  a  standard  test  problem  [6] .  This  structure  clearly  illustrates 
the  weakness  of  the  TSHELL  element.  With  a  12x12  mesh  the  error  ir 
displacement  is  17.8  percent.  The  reason  for  this  very  slov 
convergence  is  that  the  triangular  elements  essentially  add  ril 
reinforcement   to   the   very   flexible   spherical   surface. 


/• 


*/f«2.0 

(onquodront) 


fSU 


RADUS  « 10.0 
THICKNESS  "0.04 
E»  6.825  ElO7 
POISSION  RATIO  »0.3 
MESH>NxN  (onquodront) 


F«2.0 
(onquodront) 


Figure  9.   SPHERICAL  SHELL  EXAMPLE 


Table  3.   Displacement  of  Spherical  Shell 


Element  Type  and  Mesh 


Reference 

Value 

SHELL 

N  =  4 

N=8 

N=12 

TSHELL 

N=8 

N=12 

VALUE 
0.0925  in. 


0.051265 
0.075974 


ERROR 

0.0  % 

-46.1  % 

-3.3  % 

-0.7  % 

-44.6  % 
-17.8  % 
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FINAL  REMARKS 

A  new  anisotropic  triangular  shell  element,  with  normal  rotations  at 
the  nodes  has  been  developed.  The  element  can  be  connected  directly 
to   nodes   with   beam   elements   without   special   consideration. 

The  general  accuracy  of  the  element  has  been  demonstrated.  Care 
should  be  taken  if  the  element  is  used  to  model  shells  which  have 
double   curvature. 
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APPENDIX  A  -  FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 


c SHELLT 

SUBROUTINE  SHELLT (S,T,F,D,C,V,X,Y, AREA , PRES , TEMP ) 
IMPLICIT  REAL*8  (A-H,0-Z) 

c INFORMATION  CALCULATED  

C      S   =  18  X  18  STIFFNESS  MATRIX  IN  X-Y-Z  SYSTEM 

C      T   =18x18  FORCE-TRANSFORMATION  MATRIX 

C      F   =  18  x  2   FORCES  DUE  TO  PRESSURE  AND  TEMP. 

C INFORMATION  SPECIFIED  

C      D   =  6  x  6  MATERIAL  PROPERTY  MATRIX 

C  LOCAL  FORCES  ARE  ASSUMMED  TO  BE  IN  THE 

C  ORDER   Mil,  M22,  M12,  Nil,  N22  and  N12 

C      C   =  6  x  1  MATRIX  OF  THERMAL  EXPANSION  TERMS 

C      V   =  THE  DIRECTION  COSINE  ARRAY  WHICH  RELATES  THE 

C  LOCAL  1-2-3  TO  THE  GLOBAL  X-Y-Z  SYSTEM 

C      X,  Y,  Z  GLOBAL  NODE  COORDINATES 

C      AREA  =  ELEMENT  AREA 

C      PRES  =  AVERAGE  SURFACE  PRESSURE 

C      TEMP  =  AVERAGE  TEMPERATURE  CHANGE 


C 

c- 


c 


20 
C 


WRITTEN  BY  EDWARD  L.  WILSON,  JAN.  -  JUNE   1987 

DIMENSION  V (4 ,  4) , KL (24 , 3) ,XY (3) , E (6) ,H (6 , 6) , 

S(18,18),T(18,18),D(6,6),X(6),Y(6),AP(10,12) 
AM (10, 12) , A (20, 18) ,AIN(3,3) ,G (20 , 20) , JP (9) , 
JM(9) ,TT (20,18) ,F (18,2) ,C (6) 
DATA 

3,9,15/, 

12,18/ 


JP 

JM 

DATA 

1,1 

1,3 

1,2 

NT  = 


/4,5 
/1,2 
KL  / 
1,  2 
4,  7 
3,  1 
24 
CHANGE  TO 
XC  =  (  X(l) 
YC  =  (  Y(l) 
DO  20  1=1,3 
X(I)  =  X(I) 
Y(I)  =  Y(I) 
DO  20  J=l,3 
AIN(I, J) 
COMPUTE 
RED  =  0. 
AIN(1,1) 
AINU,  2) 
AIN(2,3) 
AIN(3,3) 
AINU,  2) 


10,11,16,17 
7,8,13,14,6 


2, 

10 
3, 


4 

11 

1 


4,4, 

13,14 

2,3, 


5,5,5 

17,19 
1,2,3 


20 


6 
12 

1 


6,6, 
14,1 
2,3, 


6,6,6, 

5,16,18, 

1,2,3/ 


LOCAL  COORDINATES 


X(2) 
Y(2) 


X(3) 
Y(3 


) 


)  )  / 


SYSTEM 
/  3. 
3. 


-  XC 

-  YC 


-  0.0 
INTEGRALS 
50* AREA 
=  AREA 
=-RED* ( 
=  RED* ( 
=-RED* ( 
=  AIN(2 


X(1)*X(2) 
X(1)*Y(1) 
Y(1)*Y(2) 
3) 


X(2)*X(3) 
X(2)*Y(2) 
Y(2)*Y(3) 


X(3)*X(1) 
X(3)*Y(3) 
Y(3)*Y(1) 


6. 

12. 

6. 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C FORM  "H"  ARRAY  FOR  GENERAL  6-NODE  TRIANGLE 

CALL  FORMH(H,A,X,Y) 
C FORM  "AP"  ARRAY  FOR  PLATE  ELEMENT  

CALL  FORMAP(AP,H,X,Y) 
C FORM  "AM"  ARRAY  FOR  PLANE  ELEMENT  

CALL  FORMAM(AM,H,X,Y) 
C FORM  20x20  "G"  ARRAY  

DO  160  K=l,20 

DO  160  L=1,K 
160  G(K,L)  =  0.0 
C 

DO  200  N=1,NT 

I  =  KL(N,1) 
K   =  KL(N,2) 

II  =  KL(N,3) 
DO  190  M=1,NT 
J   =  KL(M,1) 

IF  (D(I, J) .EQ.0.0)  GO  TO  190 
L   =  KL(M,2) 
IF  (L.GT.K)  GO  TO  190 
JJ  =  KL(M,3) 

G(K,L)  -  G(K,L)  +  D(I, J) *AIN(II,JJ) 
190   CONTINUE 
200  CONTINUE 
C 

DO  210  K=l,20 
DO  210  L=1,K 
210  G(L,K)  =  G(K,L) 

C FORM   20x18  "A"  ARRAY  

DO  300  J=l,18 


F(J,1)  =  0.0 

F(J,2)  =  0.0 

DO  300  1=1,20 

300 

A(I,J)  =  0.0 

DO  305  J=l,9 

J J  =  JP(J) 

DO  305  1=1,10 

305 

A(I,JJ)  =  AP(I,J) 

DO  310  J=l,9 

JJ  =  JM(J) 

DO  310  1=1,10 

310 

A(I+10,JJ)  -  AM(I,J) 



ROTATE  TO  GLOBAL  X, 

DO  350  N=l,6 

NZ  =  3*N 

NY  =  NZ  -  1 

C ROTATE  TO  GLOBAL  X,  Y,  Z  SYSTEM 


NX  =  NY  -  1 
DO  350  1=1,20 

XX  =  A(I,NX)*V(1,1)  +  A (I, NY) *V(1,2)  +  A (I ,NZ) *V (1 , 3) 

YY  =  A(I,NX) *V(2,1)  +  A(I,NY) *V(2,2)  +  A (I ,NZ) *V(2 , 3) 

ZZ  =  A(I,NX)*V(3,1)  +  A(I,NY)*V(3,2)  +  A (I ,NZ) *V (3 , 3) 
A (I, NX)  =  XX 
A  (I,  NY)  =  YY 
350  A(I,NZ)  =  ZZ 


FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 


C FORM  18x18  ELEMENT  STIFFNESS  

DO  380  1=1,20 

DO  380  J=l,18 

CALL  DOTP(G(l,I) ,A(1, J) ,SUM,20) 
380  TT(I,J)  =  SUM 
C 

DO  400  1=1,18 

DO  400  J=1,I 

CALL  DOTP(A(l,I) ,TT(1, J) , SUM, 20) 

S(I,J)  =  SUM 
400  S(J,I)  =  SUM 
C FORM  FORCE-DISPLACEMENT  TRANSFORMATION  ARRAY 

NO  =  0 

DO  500  N=l,3 

XY(1)  =  1.0 

XY(2)  =  X(N) 

XY(3)  =  Y(N) 
C 

DO  450  K=l,18 
C 

DO  410  1=1,6 
410  E(I)  =  0.0 
C 

DO  420  M=1,NT 

I  =  KL(M,1) 
L  =  KL(M,2) 
J  =  KL(M,3) 

420  E(I)  =  E(I)  +  XY(J)*A(L,K) 
C 

DO  440  1=1,6 

SUM  =0.0 

DO  430  J=l,6 
430  SUM  =  SUM  +  D(J,I)*E(J) 

II  =  NO  +  I 
440  T(II,K)  =  SUM 

C 

450  CONTINUE 
C 

500  NO  =  NO  +  6 

C FORM  THERMAL  FORCES  

IF  (TEMP. EQ. 0.0)  GO  TO  625 
DO  600  1=1,6 

CALL  DOTP(D(l,I) ,C(1) ,E(I) ,6) 
600  CONTINUE 

DO  610  1=1,6 
610  C(I)  ■  -  TEMP*E(I) 
C 

DO  620  1=1,18 
F(I,2)  =  C(1)*A(1,I) 
F(I,2)  =  C(2)*A(7,I) 
F(I,2)  =  C(3)*(  A(2,I)  +  A(6,I)  ) 
F(I,2)  =  C(4)*A(llfI) 
F(I,2)  =  C(5)*A(17,I) 
620   F(I,2)  =  C(6)*(  A(12,I)  +  A(16,I)  ) 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 


C CALCULATE  PRESSURE  FORCES 

625   IF(PRES.EQ.O.O)  GO  TO  800 

FORCE  =  -  AREA*PRES/3.0 

DO  630  1=1,18,6 

F(I   ,  1)  =  FORCE*V(l,3) 

F(I+1,1)  =  FORCE*V(2,3) 
630   F(I+2,1)  =  FORCE*V(3,3) 


800   RETURN 
c 

END 
c FORMH 

SUBROUTINE  FORMH (H, B ,X,Y) 
IMPLICIT  REAL*8  (A-H,0-Z) 
DIMENSION  H(6,6),B(6,6),Y(4),X(4) 

C FORM  COEFFICIENT  MATRIX  FOR  6  NODE  TRIANGLE  

X(4)  =  (  X(l)  +  X(2)  )  /  2. 
X(5)  =  (  X(2)  +  X(3)  )  /  2. 
X(6)  =  (  X(3)  +  X(l)  )  /  2. 
Y(4)  =  (  Y(l)  +  Y(2)  )  /  2. 
Y(5)  =  (  Y(2)  +  Y(3)  )  /  2. 
Y(6)  =  (  Y(3)  +  Y(l)  )  /  2. 


DO  100 
H(1,I) 


1  =  1 
=  1 


H(2 
H(3 

H(4 
H(5 


I) 
I) 
I) 
I) 


6 

0 

X(I) 
Y(I) 

X(I)*X(I) 
X(I)*Y(I) 
Y(I)*Y(I) 


/  2 


/  2 


100  H(6,I) 
C INVERT  TO  FORM  H  MATRIX  - 

DO  200  1=1,6 

DO  200  J=1,I 

SUM  =0.0 

DO  190  K=l,6 
190  SUM  =  SUM  +  H(I,K)*H(J,K) 

B(J,I)  =  SUM 
200  B(I,J)  =  SUM 


C 
C 


CALL  SYMSOL(B,H,6,6,0) 

RETURN 
END 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

c FORMAP 

SUBROUTINE  FORMAP (AP , H , X , Y) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AP(10,12) ,H(6,6) ,Y(4) ,X(4) ,IT(4) 

DATA  IT  /1,3,5,1/ 
C FORM  12  DOF  MATRIX  

DO  240  1=2,6 

K  =  0 

DO  240  J=l,6 

K  =  K  +  1 

AP(I-1,K)  =  0.0 

AP(I+4,K)  =  H(I,J) 

K  =  K  +  1 

AP(I-1,K)  =  -  H(I,J) 
240   AP(I+4,K)  =  0.0 
C ELIMINATION  OF  4,5,6  MID-SIDE  ROTATIONS  

X(4)  =  X(l) 

Y(4)  =  Y(l) 

DO  300  N=l,3 

DX  =  X(N+1)  -  X(N) 

DY  =  Y(N+1)  -  Y(N) 

XL  =  DSQRT(  DX*DX  +  DY*DY  ) 

S  =  DY  /  XL 

C  ■  DX  /  XL 

Tl  =  C*C/2.  -  S*S/4. 

T2  =  0.75*S*C 

T3  =  S*S/2.  -  C*C/4. 

XX  =  1.5/XL 

TC  =  XX*C 

TS  =  XX*S 
C 

11  =  IT(N) 

12  =  II  +  1 
Jl  =  IT(N+1) 
J2  =  Jl  +  1 

L2  =  6  +  N  +  N 

LI  =  L2  -  1 

I  =  N  +  6 

J  =  I  +  1 

IF  (N.EQ.3)  J  =  7 
C 

DO  300  K=l,10 

AP(K,I1)  =  AP(K,I1)  +  AP(K,L1)*T1  +  AP(K,L2)*T2 

AP(K,I2)  =  AP(K,I2)  +  AP(K,L1)*T2  +  AP(K,L2)*T3 

AP(K,J1)  =  AP(K,J1)  +  AP(K,L1)*T1  +  AP(K,L2)*T2 

AP(K,J2)  =  AP(K,J2)  +  AP(K,L1)*T2  +  AP(K,L2)*T3 

W  =  -  AP(K,L1)*TS  +  AP(K,L2)*TC 

AP(K,L1)  =  0.0 

AP(K,L2)  =  0.0 

AP(K,I)  =  AP(K,I)  +  W 

AP(K,J)  =  AP(K,J)  -  W 
300  CONTINUE 
C 

RETURN 

END 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

c FORMAM 

SUBROUTINE  FORMAM (AM, H , X, Y) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AM(10,12),H(6,6),Y(4),X(4),IT(4) 

DATA  IT  /1,3,5,1/ 
C FORM  12  DOF  MATRIX  

DO  240  1=2,6 

K  =  0 

DO  240  J=l,6 

K  =  K  +  1 


AM(I-1,K)  =  H(I,J) 

AM(I+4,K)  =  0.0 

K  =  K  +  1 

AM(I-1,K)  =  0.0 

240 

AM(I+4,K)  =  H(I,J) 

c 

ROTATE  NODE  4,  5,  6  MID-SIDE 

X(4)  =  X(l) 

Y(4)  =  Y(l) 

DISPLACEMENTS  


DO  300  N=l,3 

DX  =  X(N+1)  -  X(N) 

DY  =  Y(N+1)  -  Y(N) 
C 

NY  =  2*N  +  6 

NX  =  NY  -  1 

IX  =  IT(N) 

IY  =  IX  +  1 

JX  =  IT(N+1) 

JY  =  JX  +  1 

I  =  N  +  6 

J  =  I  +  1 

IF  (N.EQ.3)  J  =  7 
C ELIMINATE  MID-SIDE  DISPLACEMENTS  

DO  260  K=l,10 

TT  =  .125* (  -  AM(K,NX)*DY  +  AM(K,NY)*DX  ) 

AM(K,IX)  =  AM(K,IX)  +  AM(K,NX)/2. 

AM(K,IY)  =  AM(K,IY)  +  AM(K,NY)/2. 

AM(K,JX)  =  AM(K,JX)  +  AM(K,NX)/2. 

AM(K,JY)  =  AM(K,JY)  +  AM(K,NY)/2. 

AM(K,NX)  =  0.0 

AM(K,NY)  =  0.0 

AM(K,I)    =  AM(K,I)    +TT 

AM(K,J)    =  AM(K,J)    -  TT 
260  CONTINUE 
C 

300  CONTINUE 


C 


RETURN 
END 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

c LOCALT 

SUBROUTINE  LOCALT (XYZ , X , Y , V , AREA , IAXIS ) 
IMPLICIT  REAL*8  (A-H,0-Z) 
DIMENSION  XYZ(3,3) ,X(6) ,Y(6) ,V(4,4) 

C TRANSFORM  TO  LOCAL  COORDINATE  SYSTEM  

C       WRITE  (*,3000)  XYZ 
C  3000  FORMAT  (3F15.4) 

IF  ( IAXIS. NE.O)  THEN 
DO  100  1=1,3 
100      V(I,4)  =  0.0 

I  =  IABS (IAXIS) 
V(I,4)  =  IAXIS/I 
END  IF 

C DEFINE  LOCAL  1,2,3  REFERANCE  SYSTEM 

CALL  VECTOR (V (1,1) , XYZ (1,1) , XYZ (1,2) ,XYZ(1,3) , 

XYZ (2,1), XYZ (2,2), XYZ (2,3)) 
CALL  VECTOR (V (1,2) , XYZ (1,1) ,XYZ(1,2) ,XYZ(1,3) , 

XYZ (3,1), XYZ (3,2), XYZ (3,3)) 
CALL  CROSS (V( 1,1) ,V(1,2) ,V(1,3) ) 

IF  (IAXIS. NE.O)  CALL  CROSS(V(l,4) ,V(1,3) ,V(1,1) ) 
CALL  CROSS(V(l,3) ,V(1,1) ,V(1,2) ) 
C 

DO  5  N=l,3 

X(N)  =  XYZ(N,1) *V(1,1)  +  XYZ(N,2) *V(2,1)  +  XYZ (N, 3) *V (3 , 1) 
5   Y(N)  =  XYZ(N,1) *V(1,2)  +  XYZ (N, 2 ) *V ( 2 , 2)  +  XYZ  (N,  3)  *V  (3  ,  2) 

C CALCULATE  AREA  OF  ELEMENT  

AREA  =  (  X(2)*Y(3)  -  X(3)*Y(2)  +  X(3)*Y(1) 

-  X(1)*Y(3)  +  X(1)*Y(2)  -  X(2)*Y(1)  )  /  2.0 
IF  ( AREA. LE. 0.0)  THEN 
PAUSE  '  ZERO  OR  NEGATIVE  AREA  ' 
RETURN 

END  IF 
c 

RETURN 

END 
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